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ABSTRACT 

We present the Boltzmann moment equation approach for the dynamics of stars 
(BEADS-2D), which is a finite-difference Eulerian numerical code designed for the 
modelling of anisotropic and non-axisymmetric flat stellar disks. The BEADS-2D code 
solves the Boltzmann moment equations up to second order in the thin-disk approxi- 
mation. This allows us to obtain the anisotropy of the velocity ellipsoid and the vertex 
deviation in the plane of the disk. 

We study the time-dependent evolution of exponential stellar disks in the lin- 
ear regime and beyond. The disks are initially characterized by different values of 
the Toomre parameter Qs and are embedded in a dark matter halo, yielding a rota- 
tion curve composed of a rigid central part and a flat outer region. Starting from a 
near equilibrium state, several unstable modes develop in the disk. In the early linear 
phase, the very centre and the large scales are characterized by growing one-armed 
and bisymmetric positive density perturbations, respectively. This is in agreement 
with expectations from the swing amplification mechanism of short-wavelength trail- 
ing disturbances, propagating through the disk centre. In the late linear phase, the 
overall appearance is dominated by a two-armed spiral structure localized within the 
outer Lindblad resonance (OLR). During the non- linear evolutionary phase, radial 
mass redistribution due to the gravitational torques of spiral arms produces an out- 
flow of mass, which forms a ring at the OLR, and an inflow of mass, which forms 
a transient central bar. This process of mass redistribution is self-regulatory and it 
terminates when spiral arms diminish due to a shortage of matter. Finally, a compact 
central disk and a diffuse ring at the OLR are formed. An increase in Qs stabilizes the 
disks at Qs ~ 3.15, in agreement with the theoretical predictions. 

Considerable vertex deviations are found in regions with strongly perturbed mass 
distributions, i.e. near the spiral arms. The vertex deviations are especially large at 
the convex edge of the spiral arms, whereas they are small at the concave edge. The 
mean vertex deviations correlate well with the global Fourier amplitudes, reaching 
mean values of about 12° in the saturation stage. Local values of the vertex deviation 
can reach up to almost 90°. Near the convex edge of the spiral arms, the ratio of radial 
to azimuthal components of the velocity ellipsoid can deviate considerably from the 
values predicted from the epicycle approximation. 

Key words: Galaxies: general - galaxies: evolution - galaxies: spiral - galaxies: 
kinematics and dynamics 



1 INTRODUCTION 

The majority of normal disk galaxies are characterized 
by non-axisymmetric structures like spirals or bars. These 
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structural elements have been widely discussed in the lit- 
erature as a result of gravitational instabilities which are 
connected to growing density waves or global instabilities of 
disks (e.g. Binney & Tremaine 1987, afterwards BT87). 

A first insight into the properties of galactic discs was 
provided by linear stability analysis. Structure growth in a 
stellar disk within the linear regime was studied via the lin- 
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earised coUisionless Boltzmann equation by Toomre (1964), 
Zang (1976), Pichon & Cannon (1997), Evans & Read 
(1998), Jalali & Hunter (2005), and others. Toomre (1964) 
has deduced a simple stability criterion from the Jeans equa- 
tions in case of stellar disks: a stellar disk is stable against 
axisymmetric perturbations if Qs = crrrK/(3.36GE) exceeds 
a critical value Qs.crit = 1 (frr is the radial velocity dis- 
persion, K is the epicycle frequency, G is the gravitational 
constant, and E is the surface density of the disk). In case of 
non-axisymmctric perturbations the analysis becomes tech- 
nically more difficult, but the Toomre parameter turned out 
to still be a good indicator for stability, if one increases its 
critical value (Berlin et al. 1989, Polyachcnko ct al. 1997). 

A disadvantage of linear stability analysis is its restric- 
tion to small perturbations, both in amplitude and wave- 
length. Once the perturbations exceed a critical level or 
when the initial perturbations are already large, numerical 
simulations arc necessary. Hydrodynamical simulations re- 
vealed three qualitatively distinct phases for the growth of 
unstable disks (Laughlin et al. 1997). Starting from small 
amplitudes, perturbations grow first linearly. When they 
reach a critical level, mode-mode coupling becomes impor- 
tant and the non-linear evolution starts. The sclf-intcraction 
of the dominant mode might result in a rearrangement of 
the radial mass distribution (m = mode). The disk it- 
self might become violently unstable, which is indicated by 
the perturbations reaching a "macroscopic" level. Finally, 
a quasi-cquilibrium saturation stage can be reached when 
the growth of the instabilities is stopped due to shocks or 
feedback mechanisms heating up the disk dynamically. Oth- 
erwise, fragmentation continues until the disk is destroyed 
or transformed. 

Though N-body simulations have been proven to be ex- 
tremely useful for many studies of galactic dynamics, they 
are less suitable for studying the growth of very weak per- 
turbations. The main reason is that N-body calculations of 
galaxies arc usually performed with less particles than the 
actual number of stars in these systems. This causes an in- 
herent artificial particle noise which cannot be neglected for 
typical particle numbers. 

Alternatively, numerical hydrodynamics simulations be- 
came a primary tool for the analysis of growing spiral in- 
stabilities (e.g. Korchagin et al. 2000, Orlova et al. 2002). 
First, they allow for a direct comparison with the results 
of a semi-analytic linear stability analysis, including pre- 
dictions for the growth rates (Laughlin et al. 1997; 1998, 
Korchagin et al. 2000), because the same set of basic equa- 
tions is applied. Second, the standard grid methods used to 
solve the hydrodynamical equations allow to start with very 
small perturbations (basically limited by machine accuracy). 
Their evolution can be followed over many e-folding times 
until they reach the non-lineax regime and beyond. 

A major drawback of purely hydrodynamical simula- 
tions is that the dominant disk component is the stellar disk. 
The latter is characterized by an anisotropic velocity disper- 
sion, whereas the hydrodynamical equations deal with an 
isotropic pressure. Secondly, one has to assume an equation 
of state connecting pressure, (surface) density and temper- 
ature (velocity dispersion). For instance, if one assumes a 
polytropic-like equation of state, i.e. P = CT,'', the set of 
hydrodynamic equations is closed with the Euler equation 
describing the momentum transport (for an example of a star 



bility analysis based on such assumptions see e.g. Aoki et al. 
1979). Since the temperature does not show up in the poly- 
tropic equation of state, no energy transport equation is re- 
quired. Obviously, this is a convenient simplification which, 
however, might not be generally applicable. Moreover, the 
profile of the Toomre parameter is completely fixed by the 
surface density distribution. 

A more complete ansatz for describing stellar disks is 
the coUisionless Boltzmann equation (BT87) 



^-0 
dt dr dv 



(1) 



where f{r,v,t) is the distribution function. The term v in- 
cludes accelerations like gravitational forces, pressure terms, 
etc. Since a direct solution is already numerically very dif- 
ficult, the problem is often reduced by taking the velocity 
moments of equation (1). By this, one gets an infinite se- 
ries of moment equations where the equations for moment 
i include moments of the order i -|- 1. The knowledge of all 
velocity moments is equal to the knowledge of the distri- 
bution function itself. Practically, one has to terminate the 
set of equations by a closure condition. Setting all moments 
of third order to zero (zero-heat-flux condition) results in 
the Jeans equations which are well applicable to stellar sys- 
tems with negligible two-body relaxation. In that case, the 
second order moment equations which yield information on 
the velocity ellipsoid are included and stellar anisotropy can 
be considered. An example are the one- and two-dimensional 
chemo-dynamical models by Theis et al. (1992) and Samland 
et al. (1997). For the special case of galactic disks, Amendt 
& Cuddoford (1991) and Cuddeford & Amendt (1991) stud- 
ied the higher order moments in detail. They found that 
for reasonable constraints to the distribution function the 
third order terms vanish to leading order in the plane of the 
disk. This corroborates the assumption of a zero-heat-flux. 
Extensions of the Boltzmann moment equations including 
coUisional processes have also been used (e.g. Larson 1970, 
Louis 1990 or Giersz & Spurzem 1994 for the evolution of 
star clusters). 

In this paper we present the Boltzmann equation ap- 
proach for the dynamics of stars (BEADS-2D), which is a nu- 
merical code that solves the Boltzmann moment equations 
up to second order using the methods of finite-differences. 
The BEADS-2D code is applied to study the time-dependent 
evolution of flat, non- isotropic stellar disks in the linear 
regime and beyond. In contrast to previous numerical stud- 
ies, we include the non- axisymmetric velocity dispersion 
terms into the Boltzmann moment equations. Therefore, our 
approach allows us to study the evolution of the anisotropy 
and the vertex deviation of the stellar component. We com- 
pare the stability properties of non-isotropic stellar disks 
with the predictions of linear stability analysis for non- 
isotropic disks by Polyachenko et al. (1997). We provide a 
physical interpretation for growing instabilities in our model 
stellar disks. 

In the following section we describe the Boltzmann 
moment equation approach. The basic principles of the 
BEADS-2D code are discussed in § 3. The results of nu- 
merical modelling are presented in § 4. Finally, conclusions 
and a summary are given in § 5. 
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2 STELLAR SYSTEMS: THE MOMENT 
EQUATION APPROACH 

2.1 The 3D-Boltzmann moment equations 

Stars are collisionless objects that move on orbits deter- 
mined by the large-scale gravitational potential $. An ex- 
act description of such a system requires the solution of 
the collisionless Boltzmann equation (1) for the distribu- 
tion function of stars f{r,v,t) in the phase space {r,v). 
The distribution function of stars / = f{r,v,t) contains 
a fundamental description of the stellar system. Its lower 
order moments like the density p = f fd^v, the mean ve- 
locity u = j f V dPv, or the velocity-dispersion tensor 
(^Ij = / / ~ — Uj) dPv of stars can be deduced 

- at least partly - from observations. The corresponding 
equations for p, pu, and pafj can be obtained by taking 
moments of the collisionless Boltzmann equation (1). If one 
closes the system by assuming the zero-heat-flux approxi- 
mation Qijk = p~^ J f {vi— Ui){vj — Uj){vk - Uk) (fv = 0, a 
set of ten partial differential equations (the velocity disper- 
sion tensor is symmetric by definition) can be derived (see 
Appendix A). 

The three off-diagonal elements of the velocity- 
dispersion tensor can be regarded as a measure of alignment 
of the principal axes of the velocity ellipsoid with respect to 
the coordinate system. Wc note that these terms may bo- 
come both positive or negative (although they have a square 
mark) , as follows directly from the definition of the velocity- 
dispersion tensor. 

In our models, the gravitational potential $ is composed 
of the sum of an external contribution $cxt (r, t) and the self- 
gravity 3>sys(r,f) of the disk mass distribution p(r) obeying 
the moment equations. Its potential can be derived from the 
Poisson equation 



A$syB = 47rGp. 



(2) 



The external term describes amy known (eventually time- 
dependent) potential like a dark matter halo, whereas the 
"system" corresponds to the live disk. 



Continuity equation: 



dt 



i|-(r.S.«.) + iA(E.. 
r or r 0(p 



= 0. 



(3) 



Momentum equations: 
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Velocity dispersion equations 
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2.2 The 2D-Boltzmann moment equations for flat 
disks 

A substantial fraction of stars in spiral galaxies are concen- 
trated in the disk. Hence to a first approximation, stellar 
disks may be regarded as having zero thickness. In the 
thin-disk approximation, all motions are localized within 

the (r, (f)) plane. The moment equations in the thin-disk 
approximation can be obtained from equations (Al)-(AIO) 
by assuming a vanishing z-gradient of all physical variables 



and setting Uz = (negligible vertical motion), 



0, 



and uIj, = 0. The stellar volume density p is vertically 
integrated to yield the surface density E. The assumptions 
of air = 0, and cr^^ = are justified on observational 
grounds (Binney & Merrifield 1998). The vertical velocity 
dispersion Ozz is obtained by assuming a constant ratio of 
the vertical to the radial (azimuthal) velocity dispersions. 
Below, wc provide the moment equations in the thin-disk 
approximation written in polar coordinates. For the conve- 
nience of coding, the advection terms are set in brackets. 



^(EfT^,^) -I- 



15/ ^2 . 1 d 2 



1 du4, 1 
r d(f> j 
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= 0. 



(8) 



According to our experience, equations (3)-(8) are expressed 
in the computationally most stable form. For convenience we 
will use the name kinetic equations or kinetic models when 
we solve the second order Boltzmann moment equations. 



2.3 The gravitational potential 

In our models the gravitational potential consists of two 
parts, a live disk and a static halo component. The grav- 
itational potential $disk of the flat disk can be calculated by 
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(BT87, Sect. 2.8) 



$diBk(?-, ( 



poo 

= -G r'dr' 
Jo 

^ p E(r',0')#' 
Jo 



^ + J.2 _ 2rr' cos(0' — ^) 



(9) 



This sum is calculated using a FFT technique which applies 
the 2D Fourier convolution theorem for polar coordinates. 

The halo properties are fixed by the rotation curve 
parametrized by 



— ) 

rflat/ 



(10) 



1 + 



The transition radius between an inner region of rigid ro- 
tation and a flat rotation in the outer part is given by rsat 
which we set to 3 kpc. The smoothness of the transition is 
controlled by the parameter n^, set to 2. The velocity at 
infinity, Woo, was set to 208 km s~^. The solid line in Fig. 1 
shows the rotation curve we used. The corresponding halo 
potential $haio is then derived from equation (4), where 3> 
is substituted with $haio + $disk and Ur is set to zero. 



3 NUMERICAL MODELLING 
3.1 The BEADS-2D code 

The 2D-equations for thin disks arc discrotizcd on an Eulc- 
rian grid with a logarithmic grid spacing in the r-direction 
and uniform grid spacing in the (^-direction. The differ- 
ent terms are taken into account by applying an operator 
splitting technique similar to the ZEUS program (Stone & 
Norman 1992). Advcction is performed by the second-order 
van-Leer scheme. The time step is determined according to 
the usual Courant-Priedrichs-Levy criterion modified to take 
into account the velocity dispersion. Namely, we diagonaJize 
the velocity dispersion tensor at each time step and use the 
diagonal components in the time-step limiter. For our simu- 
lations we use reflecting boundary conditions in the radial di- 
rection (and periodic boundary conditions in the azimuthal 
direction). Our simulations are done on a grid with 256 x 256 
cells covering the radial range from 0.2 kpc to 30-45 kpc (de- 
pending on the model). The gravitational force of a thin disk 
near its inner boundary at r = 0.2 kpc is directed radially 
outward, which is an artificial effect due to the lack of ma- 
terial within the inner boundary (note that such a problem 
would not exist in the case of spherical symmetry) . In order 
to reduce the magnitude of this spurious outward gravity 
force, we introduce an inner circular disk of constant den- 
sity and outer radius Tout = 0.2 kpc. This inner disk merely 
serves as a central gravity source and its gravitational po- 
tential $id can be computed by decomposing the inner disk 
into a series of concentric circular rings of constant density 
(see Vorobyov & Basu 2006 for details). A more accurate 
method that is used in this paper involves an expansion of 
$id into Lcgcndrc polynomials in the plane of the disk. The 
details of this method will be given in a future paper. The 
total gravitational potential $ = $disk + *haio + ^id is then 
used in the momentum equations. 



3.2 Artificial viscosity 

The artificial viscosity, much like kinematic viscosity in a 
real fiuid, is often used to smooth discontinuities where the 
finite-difference equations break down. The 2nd order Boltz- 
mann moment equations have to be modified to account for 
the viscous stresses and dissipation due to the artificial vis- 
cosity. 

The implementation of the artificial viscosity in the ki- 
netic model requires the use of a tensor formalism. The ar- 
tificial viscosity stress tensor Q can be written by analogy 
to the general form of the molecular viscosity stress tensor 



Q = 2A*v[Vw--(V-M)e], 



(11) 



where Vtt is the symmetrized velocity gradient tensor, V • u 
is the velocity divergence, e is the unit tensor, and /Ltv is 
coefficient of artificial viscosity defined as 







if V • u < 
otherwise. 



where L is a constant with dimensions of length, chosen to 
be one zone width. The components of the artificial viscosity 
stress tensor Q in cylindrical coordinates are given in Ap- 
pendix B. The off-diagonal components of Q are often set 
to zero to ensure that the artificial viscosity smoothes only 
compressible shock fronts and leaves large gradients of shear 
unchanged. The corresponding viscous stress and dissipation 
terms are introduced into the kinetic equations (4)-(8) by 
defining a generalized stress tensor Pij = T,aij -\- Qij , where 
the artificial viscosity stress tensor Qij is analogous to the 
viscous stress tensor —nij in gas dynamics. We note that 
Qij has a plus sign because the coefficient of artificial vis- 
cosity /iv is negative by definition. The modified momentum 
equations can be obtained from equations (4) and (5) by 
substituting "Earr with Prr and Ecrg,^ with P^tp. The modi- 
fied velocity dispersion equations can be derived by applying 
the same substitution to the last two terms of equations (6) 
and (7) and to the last three terms in equation (8). 



3.3 Units 

The units are chosen to be 10^" Mq and 1 kpc for the mass 
and length scales, respectively. The gravitational constant G 
was set to 1. This results in a velocity unit of 207.4 kms"'^ 
and a unit for the angular speed of 207.4 kms~^ kpc~^. The 
time unit is 4.7 Myr. These units are used throughout the 
paper, unless other units are given explicitly. 



3.4 Tests and accuracy 

We performed several tests for our numerical code. Two 
standard tests - one focusing on the accuracy of the im- 
plementation of the advection and another one on the con- 
servation of the specific angular momentum - are given in 
Appendix C in some detail. Both tests demonstrate the abil- 
ity of the code to handle the advection terms numerically 
well, so that we do not expect spurious structure formation 
due to advection problems. 
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Figure 1. Initial profiles of the rotation curve (solid line, identical 
for all models), and Toomre parameter (dashed line, model K3) 
We note that the Q-profile just scales in diflferent models. 



Table 1. Model parameters of stellar disks 



Model 


Qs 


CTrr {r = 0) 


ilout (kpc) 


Kl 


1.1 


112.5 km s-l 


30 


K2 


1.3 


133 km s-l 


30 


K3 


1.6 


164 km s-i 


45 


K4 


2.5 


250 km s-i 


45 


K5 


3.15 


315 km s-i 


45 



3.5 Fourier amplitudes 

In order to visualize the nature and growth of instabilities 
in disks, it is common to use the azimuthal Fourier modes 
Om at a given radius r. They are calculated for m ^ by 
(e.g. LaughUn et al. 1998) 



am(r,t) = ^ 



T,{r,cf>,t)e'""^d4> 



(12) 



An even simpler measure for the global structure of a disk 
are the (radially integrated) global Fourier modes defined by 

C^{t) = —- / a^{r,t)rdr. (13) 

Aidisk J^^^^ 

Mdisk is the mass of the disk within the radial range 
[riow, '"high]- The global Fourier amplitude Am is the mod- 
ulus of Cm, whereas the global growth rate is given by the 
time derivative of the logarithmic Fourier amplitude, i.e. 



7" 



d(hi|C.n(0 
dt 



(14) 



4 RESULTS 

4.1 Stability properties of non-isotropic (stellar) 
disks 

The analysis of the stability of a thin stellar disk to a local 
axisymmetric perturbation states that the disk is stable if 
(Toomre 1964) 



According to Polyachenko et al. (1997), the Toomre criterion 
(15) should be modified to Qs > 27r/?/3.36, if local non- 
axisymmetric perturbations are considered. Here, /3 is a 
function of the rotation curve. In the most interesting case 
of a fiat rotation curve, /3 = 1.69, so that a thin stellar disk 
becomes stable if Qs ^ 3.15. 

In this section we numerically study the stability prop- 
erties of thin stellar disks characterized by different values of 
Qs- We do not focus on the detailed study of growth rates 
of unstable modes in the linear regime, since it would re- 
quire a detailed comparison with previous analytical lineax 
stability analyses. Instead, we simply compare the stabil- 
ity properties of stellar disks obtained using the BEADS- 
2D code with the analytical predictions of Polyachenko et 
al. (1997). By using the kinetic equations (3)-(8), we take 
full account of the inherent anisotropic properties of stellar 
disks. We present the results for five models, the parameters 
of which are shown in Table 1. The acronym K stand for 
"Kinetic models". The outer radii of stellar disks i?out were 
chosen so that to exclude the infiuence of the outer refiecting 
boundary. The surface density distribution and the rotation 
curve are identical for all models, while the initial velocity 
dispersion distributions are different. The surface density is 
distributed exponentially according to 



E(r) 



(16) 



with a radial scale length rj, of 4 kpc. The central surface 
density Eo is set to lO'"* Mq pc^^. Thus, the total mass of 
the disk within the 30 kpc radius is approximately IO^^Mq. 
We note that the densities of stellar disks are indeed found 
to decay exponentially with distance, with a characteris- 
tic scale length increasing from 2-3 kpc for the early type 
galaxies to 4-5 kpc in the late type galaxies (Freeman 1970). 
The initial azimuthal velocity is chosen according to equa- 
tion (4), whereas the radial velocity vanishes initially. The 
radial component of the velocity dispersion is obtained from 
the relation arr = 3.36QsGE/k for a given value of the 
Toomre parameter Qs. We assume that throughout most of 
the disk Qs is constant (and equal to the value given in Ta- 
ble 1) but is steeply increasing with radius at r > 30 kpc. 
The initial radial distribution of Qs in model K3 is shown in 
Fig. 1 by the dashed line. The azimuthal component of the 
velocity dispersion a^^ is determined adopting the epicycle 
approximation, in which the following relation between cr^^ 
and a^r holds (Binney & Tremaine 1987, p. 125): 



du^ 
dr 



+ 1 



(17) 



Qi 



3.36GE 



> 1. 



(15) 



The off-diagonal component is initially set to zero. In 
the beginning of simulations we introduce a small random 
perturbation with a maximum relative amplitude 10"^ to 
the initial density distribution. 

Figure 2 shows the temporal behaviour of the dominant 
m = 2 mode in models K1-K5 with different initial values 
of the Q-parameter. It is evident that models with initially 
larger Qs attain the saturation phase at later times than 
do the models with smaller Qg. For instance, in model K2 
(Qs = 1.3) the saturation of the dominant mode is reached 
after 1.4 Gyr, whereas in model Kl (Qs = 1.1) the dominant 
mode saturates at an earlier time t = 0.8 Gyr. Model K4 
with the initial Qs = 2.5 reaches the saturation phase only 
after approximately 7 Gyr. The global Fourier amplitudes in 
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Figure 2. The temporal behaviour of the m = 2 mode in kinetic 
models with different initial Q-parameters: Qs = 1.1 (Kl), 1.3 
(K2), 1.6 (K3), 2.5 (K4), and 3.15 (K5). 




10 15 20 

Radial distance (kpc) 



30 



Figure 3. The initial radial distributions of the X-paramctcr for 
the m = 1 mode (dashed line) and the m = 2 mode (solid line). 
Note that X is identical in all five models listed in Table 1, be- 
cause these models initially differ only in the values of the velocity 
dispersion. 



the saturation phase arc noticeably larger in models Kl K3 
(C2(t) ^ 0.1 - 0.25) than in models K4 and K5 (C2(i) « 
0.02 — 0.05), indicating a stronger instability in colder disks. 

The linear stability analysis of non-isotropic stellar 
disks to local perturbations performed by Polyachenko et 
al. (1997) yields a critical value of Qs.crit ~ 3.15 for stellar 
disks with a flat rotation curve. Our numerical simulations 
indicate that the m = 2 mode in the Qs = 3.15 model K5 
saturates at a very low level of C2 « 0.02 — 0.03. As we shall 
see in § 4.3, stellar disks characterized by such small values 
of C2 can only develop a clumpy pattern with very small 
positive density perturbations rather than a spiral pattern 
or a bar. 
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Figure 4. The positive density perturbation that develops in the 
inner 10 kpc region in model K2 at t = 1.0 Gyr. The lumpy struc- 
ture is naturally explained as the result of interference between 
leading and trailing short-wavelength disturbances that propa- 
gate through the centre of the disk. The scale bar is in Mq pc~^. 



4.2 Swing amplification 

An appealing physical interpretation for the growth of 
strong instabilities that are common in numerical simula- 
tions of differentially rotating stellar disks is swing amplifica- 
tion (Toomre 1981). Amplification occurs when any leading 
disturbance unwinds into a trailing one due to differential 
rotation. It is helpful to introduce the parameter X = A/Acr 
when discussing the efficiency of swing amplification, where 
A = 2Trr/m is the circumferential wavelength of an m-armed 
disturbance and Acr = 47r^GE/K^ is the longest unstable 
wavelength in a cold disk. According to Julian and Toomre 
(1966), in a stellar disk with a fiat rotation curve the gain 
of the swing amplifier is the largest when 1 < X < 3. The 
gain also strongly depends on the value of the Q-parameter. 
Substantial amplification is also possible for X < 1 if Q is 
safely below 2. 

Figure 3 shows the initial radial distribution of the X- 
parameter for the m = 1 (dashed line) and m = 2 (solid 
line) disturbances. It is evident that the best conditions for 
swing amplification (1 < X < 3) of the m — 1 mode are met 
in the inner portion of the disk at (0.85 — 4.7) kpc, while the 
swing amplifier of the m = 2 mode is maximal at intermedi- 
ate radii (4.7— 12.8) kpc. However, feedback loops that turn 
leading disturbances into trailing ones must be present in 
stellar disks in order for swing amplification to destabilize 
the disk (see e.g. BT87 for a discussion). Since our stel- 
lar disks (initially) have no inner Lindblad resonances, the 
trailing disturbances can propagate through the centre and 
emerge on the other side as leading ones, thus providing a 
feedback for the swing amplifier. Figure 4 supports this point 
of view. It shows the positive density perturbation that de- 
velops in the inner 10 kpc in model K2 after 1 Gyr from the 
beginning of simulations. The lumpy structure that is seen 



Boltzmann moment equation approach for stellar dynamics 7 



7 6 6 4 3 2 1 32 28 25 21 18 14 11 7 4 





20 


o 




Q. 






10- 












in 




c 





% 








as 


-10- 


T3 


CD 




CC 




-20^ 



1O0O Myr(-2.5) 


ftSflUdyrf-I.e) 















-la 10 20 -20 -10 

Radial distance {kpc) 



« as 



14$- t3B 114 98 82 



320 28S 250 215 180 145 110 7i 40 5 



^ 20- 

o. 



CO 

■6 

^ -lOH 



T3 

a: 



•f30aMyr(-1.1) 


1400 Myr (-0.63) 






















-10 6 io 20 


-20 --io a 10 a'a 



Radial distance (kpc) 

Figure 5. Positive density perturbations in model K2 at four 
different times indicated in each panel. The quantities in paren- 
theses give the values of C2 (t) (in log units) at the corresponding 
evolutionary times. The outer and inner dashed circles sketch the 

positions of corotation and the outer Lindblad resonances. The 
scale bars give the positive density perturbations in Mq pc~^. 



in Fig. 4 is naturally explained as the result of interference 
between leading and trailing short-wavclcngth (A < Acr) 
disturbances, propagating through the disk centre. A sim- 
ilar phenomenon was seen in the numerical modelling by 
Toomrc (1981) and was discussed by BT87 (cf. Sect. 6.3, 
2(c)). The m = 1 density response dominates the innermost 
disk, while the m = 2 density perturbation is the strongest 
in the intermediate and outer regions in Fig. 4. This ten- 
dency can indeed be expected from the radial distribution 
of the X-parameter for the m = 1 and m = 2 modes shown 
in Fig. 3. The perturbation amplitudes gradually decline at 
larger radii. 

To better illustrate the development of the m = 2 insta- 
bility, we plot in Fig. 5 the positive density perturbation ob- 
tained in model K2 at four different times. The quantities in 
parentheses give the global Fourier amplitudes (in log units) 
at the corresponding times, while the scale bars give the am- 
plitudes of positive density perturbations in Mq pc~^. It is 
clearly seen that the m = 2 mode is the fastest growing 
mode. The lumpy structure ai t = 1 Gyr is gradually trans- 
forming into an m = 2 spiral perturbation at t = 1.15 Gyr. 
The spiral perturbation is mostly localized within the outer 
Lindblad resonance (OLR; shown by the outer dashed cir- 
cle) , which is in good agreement with the theoretical predic- 
tions. 

As the global Fourier amplitude C2{t) exceeds 0.1 at 
t « 1.3 Gyr, the subsequent evolution of the m = 2 spi- 
ral disturbance is governed by the non-linear effects of mass 
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Figure 6. The aaimuthally averaged surface density distribution 
of stars (thin solid and dashed lines) and the azimuthally averaged 
gravitational torque (thick solid and dashed lines) obtained in 
model K2 at two different evolutionary times as indicated in the 
legend. 
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Figure 7. Positive density perturbations in model K2 at 1.5 Gyr 
(left panel) and 1.75 Gyr (right panel). The scale bar is in 

Mq pC-2. 



and angular momentum redistribution in the disk. Spiral dis- 
turbances are known to transfer mass inward and angular 
moment outward (see e.g. Laughlin et al. 1997). The az- 
imuthally averaged gravitational torque r(r), which is the 
sum of the individual torques r = — m9$/90 exerted on all 
computational cells at a given radius r, is known to be a 
good diagnostic tool for the mass and angular momentum 
redistribution in the disk. Figure 6 shows r(r) and the az- 
imuthally averaged surface density distribution in model K2 
obtained at t = 1.3 Gyr and t — 1.4 Gyr. The stars that are 
characterized by negative F are losing angular momentum 
and spiralling into the centre, while the stars that are dis- 
tinguished by positive F are gaining angular momentum and 
moving radially outward. As a result, the density in the cen- 
tral region starts to grow considerably after 1.3 Gyr, which 
triggers the formation of a transient bar inside corotation. 
This process is evident in the right lower panel of Fig. 5. 
The material that flows radially outward settles into a weak 
outer resonant ring near the position of the outer Lindblad 
resonance at approximately 15 kpc (note that F becomes 
zero outside the outer Lindblad resonance). 

A shortage of matter that develops in the intermediate 
region weakens the spiral structure and terminates the mass 
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Table 2. Pattern speeds and position of resonances 



Model 




CR 


OLR 


Kl 


30.2 


6.2 


11.5 


K2 


23.1 


8.4 


15 


K3 


19.1 


10.4 


18.1 


K4 


11.5 


17.6 


24.1 



the pattern speeds are in km s"'^ kpc""' and positions of 
corotations (CR) and outer Lindblad resonances (OLR) are in 
kpc. 

and angular momentum redistribution. In the end, the bar 
transforms into a compact dense central disk, while the spi- 
ral structure virtually disappears. The outer resonant ring 
occupies a substantial portion of the outer disk. Figure 7 
illustrates this process and shows the positive density per- 
turbation that develops in model K2 at t = 1.5 Gyr (left) 
and t = 1.75 Gyr (right). This phenomenon is a nice exam- 
ple of the dominant mode self-interaction originally studied 
by Laughlin et al. (1997) in the context of gaseous disks. 



4.3 Stabilization at larger Qs 

The stabilizing effect of increasing Qs is clearly seen in Fig. 2. 
The m = 2 perturbations in hotter disks are characterized 
by lower growth rates and they take longer time to satu- 
rate. The global Fourier amplitudes in the saturation phase 
are noticeably larger in models K1-K3 (C-zit) ~ 0.1 — 0.25) 
than in models K4 and K5 (C2(t) ^ 0.02 ~ 0.05), indicating 
a weaker instability for hotter disks. Figure 8 shows the pos- 
itive density perturbation that develops in model K4 in the 
saturation phase at 6.8 Gyr (left panel) and 7.0 Gyr (right 
panel). The bar and the spiral arms have a patchy, flocculent 
appearance. The positive density perturbation (in physical 
units of Mq pc~^) in the saturation phase is on average 
an order of magnitude smaller than in model K2 (Fig. 5). 
Clearly, the growth of the global modes is considerably sup- 
pressed in model K4. Figure 9 shows the positive density 
perturbation that develops in model K5 in the saturation 
phase at 7.6 Gyr (left panel) and 8.8 Gyr (right panel). It is 
obvious that the positive density perturbation in model K5 is 
characterized by a clumpy structure of low amplitude rather 
than by a regular spiral or bar. This is because in model K5 
both the m = 1 and m = 2 modes saturate at a nearly equal 
(and very small) value of Ci,2 « 0.02 — 0.03. Thus, we con- 
firm the theoretical predictions of Polyachenko et al. (1997), 
who have analytically shown that stellar disks characterized 
by Qs < 3.15 and a flat rotation curve are unstable to local 
density perturbations. 

We have shown in Sect. 4.2 that the physical interpre- 
tation for the instability growth in our models is swing am- 
plification, with the feedback provided by the trailing dis- 
turbances propagating through the disk centre. Julian and 
Toomre (1966) have shown that the gain of swing amplifi- 
cation is greatly sensitive to the values of the X- and Q- 
parameters. More specifically, the gain diminishes consid- 
erably for Qs > 2. Since the initial radial profiles of the 
X-parameter are identical in all model disks (see Fig. 3), 
the increase in Qs and, by implication, in random motions 
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Figure 8. Positive density perturbations in model K4 at 6.8 Gyr 
(left panel) and 7.0 Gyr (right panel). The quantities in paren- 
theses give the values of C'2 (t) (in log units) at the corresponding 
evolutionary times. Note a patchy, flocculent structure of the bar 
and spiral arms. The scale bar is in Mq pc~^. 
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Figure 9. The same as in Fig. 8 but for model K5. 

of stars suppresses the short-wavelength disturbances and 
thus inhibits the swing amplifier. 

The comparison of dynamical properties of our model 
disks suggests another stabilization effect that may be 
present in hotter stellar disks. Table 2 gives the pattern 
speeds and the positions of corotation and outer Lindblad 
resonances in models K1-K4. These values are computed 
in the saturation phase and may slightly change during the 
evolution. It is clearly seen that models with larger Qs are 
distinguished by a slower pattern speed f2p. This is in qual- 
itative agreement with the predictions of the linear stability 
analysis of round galactic disks by Pichon & Cannon (1997), 
who have shown that colder disks yield more centrally con- 
centrated and faster rotating spirals. Slower pattern speeds 
in hotter disks imply that the latter can easier develop the 
inner Lindblad resonance due to the rearrangement of mass 
and angular momentum. Indeed, solid lines in Fig. 10 show 
the initial radial profiles of stellar angular velocity Q and 
quantities fl ± k/2, while the dashed lines give the pattern 
speeds Op in models K1-K4 (from top to bottom). None 
of the models have the inner Lindblad resonance but the 
tendency is clear: hotter disks have better chances for the 
inner Lindblad resonance to occur. The inner Lindblad reso- 
nance would prohibit the propagation of spiral disturbances 
through the disk centre, thus providing the cutoflf for the 
swing amplification feedback. 
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Figure 10. Radial profiles of tlic stellar angular velocity Q 
(dotted-dashed line) and quantities Q±k/2. The pattern speeds 
Up of models K1-K4 are shown by the dashed lines (from top to 
bottom). 



4.4 The m = 1 density perturbation 

The initial radial distributions of the X-parameter in Fig. 3 

(identical for all our model disks) suggest that the m — 1 
perturbations should dominate in the innermost parts of 
stellar disks. In fact, this can already be seen in Fig. 4, 
where the positive density response in the inner 1 kpc has a 
crescent shape. We find that hotter disks are more strongly 
disposed to the development of the m — 1 density pertur- 
bations than colder disks. Although the instability in disks 
with Qs < 2.5 is dominated by the symmetric m = 2 mode, 
the lopsided m = 1 mode in disks with Qs > 2.5 may 
compete with or even prevail over the m = 2 mode. This 
tendency is clearly scon in Fig. 11, in which wc show the 
temporal evolution of the m = 1 mode (dashed lines) and 
m = 2 mode (solid lines) in models K3-K6. Model K6 is 
initially characterized by Qs = 3.5. It is obvious that the 
m = 2 mode dominates during the evolution (except for 
the very early phase) in models with Qs = 1.6 (K3) and 
Qs = 2.5 (K4). This tendency is not seen in the Qs = 3.15 
model K5, in which both the m = 1 and m = 2 modes have 
a nearly equal global Fourier amplitudes and growth rates 
in the saturation phase. Finally, the Qs = 3.5 model K6 is 
characterized by the dominant m = 1 mode. 

The preponderance of the m = 1 mode in model K6 is 
due to the fact that one-armed perturbations are much more 
difficult to stabilize (Evans & Read 1998). Indeed, the m — 1 
modes have no inner Lindblad resonance (the quantity Q — n 
is negative), which removes a powerful stabilizing effect for 
the m = 1 mode. On the other hand, the stabilizing effect of 
the inner Lindblad resonance for the m = 2 mode is expected 
to become stronger along the sequence of increasing Qs (sec 
Fig. 10). As a consequence, the strength of the m = 2 mode 
in the saturation phase (as determined by C2{i)) decreases 
for hotter disks, whereas the strength of the m = 1 mode 
(as determined by Ci(t)) stays nearly constant in all four 
models shown in Fig. 11. We note that according to Evans 
& Read (1998) stellar disks with a sharp central cut-out 
(stellar density drops to zero in the centre) are much more 
unstable to the development of the m = 1 modes. 
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Figure 11. Temporal evolution of the m = I modes (dashed 

lines) and the m = 2 modes (solid lines) in models K3, K4, K5, 
and K6. Model K6 is characterized by Qs = 3.5. Note that in 
the saturation phase the global Fourier amplitudes of the m = 2 
modes decrease along the sequence of increasing Qs, while the 
global Fourier amplitudes of the m = 1 modes stay nearly con- 
stant in all four models. 



4.5 Vertex deviation 

The non-isotropic nature of the BEADS-2D code allows for 
a numerical study of the velocity ellipsoid properties in spi- 
ral galax;ies. The non-eixisymmetric component of the spiral 
gravitational field may produce considerable perturbations 
on stellar orbits. A vertex deviation ly can be used to study 
the magnitude of the spiral gravitational field (Kuijken & 
Tremaine 1994). The vertex deviation has a simple geomet- 
ric interpretation: it is the angle between the major axis of 
the velocity ellipsoid (the surface that has semi-axes azz, 
(Jrr, and a^cf,) at a given position in the disk and the centre- 
anticentre direction. Based on the above explanation, the 
value of Zv can be defined as (Binney and Merrifield 1998, 
p. 630) 



Zv = -atan 



2o-^^ 



(18) 



This definition is incomplete because it implicitly assumes 
the epicycle approximation, i.e a^r > c^ltt,- As a result, the 
value of Zv is limited to the [—45°, 45°] range depending 
on the sign of cr^^. As a matter of fact, the epicycle ap- 
proximation may break down in the presence of a strong 
non-axisymmetric gravitational field of spiral arms. Con- 
sequently, the azimuthal dispersion may become (locally) 
larger than the radial one and the angle between the major 
axis of the velocity ellipsoid and the centre-anticentre direc- 
tion may exceed ±45°. To take this possibility into account, 
we extend the definition of the classical vertex deviation as 
follows: 



Zv = 




> 



4"l> 



Black dashes in Fig. 12 show the major axes of the veloc- 
ity ellipsoids superimposed on the positive density perturba- 
tion map which is obtained in model K2 at t = 1.4 Gyr. Con- 
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siderable vertex deviations are seen near the convex edges of 
spiral perturbations and near the bar, where Iv may become 
as high as 89° . Interestingly, the concave edges of spiral per- 
turbations do not produce noticeable vertex deviations. The 
outer axisymmctric regions and the intor-arm regions arc 
also characterized by near-zero vertex deviations. The tem- 
poral evolution of a mass- weighted value of ly over the entire 
disk is approximately correlated with C2{t). For instance 
Zv = 1.8° at t = 1.0 Gyr and Zv = 12.6° at t = 1.4 Gyr. 

Since wc start our simulations with a zero vortex de- 
viation {artf, = 0), we confirm that the spiral gravitational 
field produces a vertex deviation - a conclusion previously 
made by Kuijkon & Trcmaine (1994) on analytical grounds. 
A more detailed study of the properties of the velocity el- 
lipsoids in spiral galaxies will be presented in a forthcoming 
paper. 

In contrast to gaseous disks, stellar disks are essentially 
non-isotropic. In the epicycle approximation, the ratio of 
radial to azimuthal components of the velocity ellipsoid can 
be described by equation (17), implying that stellar disks 
are colder in the azimuthal direction than in the radial di- 
rection. The initial ratios cr,^^ : ffrr in model K2 lie in the 
(1 — 0.71) range, with the minimum and maximum values 
being near the centre and at the outer edge of the stellar 
disk, respectively. As the disk evolves and the spiral struc- 
ture emerges, the ratios cr^^ : (Trr start to show considerable 
deviations from the initial values. In the saturation phase, 
the bulk of the stellar disk is characterized by : arr lying 
in the 0.65 — 0.75 range, which is comparable to a measured 
mean value of a^^ : arr ~ 0.8 for the disk of NGC 488 
(Gerssen et al. 1997). However, considerably smaller values 
of (7^0 : cjrr < 0.4 are found at the leading edges of spi- 
ral arms, implying that the epicycle approximation breaks 
down there. A more detailed study of the anisotropy of the 
velocity ellipsoid will be presented in a future paper. 



4.6 Miscellaneous 

In all models considered in this paper, the outer reflecting 
boundary is set at a large enough distance (r > 30 kpc) so 
as to exclude its influence on spiral generation and growth. 
A factor of 2 variation in the position of the inner reflect- 
ing boundary (0.2 kpc in the present simulation) does not 
introduce a noticeable change in the numerical results. The 
radial and azimuthal components of velocity dispersions may 
accidentally become negative in the late saturation phase, 
probably due to large gradients of the stellar potential. In 
that case the time step is reduced and the solution is sought 
again. In the rare case that this procedure fails to resolve the 
problem, we set the velocity dispersion in particular compu- 
tational cells to a predefined small value. 

Our simulations have shown that artificial viscosity is 
usually not needed in models with Qs > 1.3. However, the 
artificial viscosity has to be included if the evolution of spiral 
instabilities is followed deeply into the nonlinear regime in 
models with Qs < 1.3. The solution is only slightly modified 
by the artificial viscosity. For instance, the relative difference 
in the growth rate of the dominant m = 2 mode in model Kl 
with and without artificial viscosity is kept below 2% during 
the simulations. 
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Figure 12. Tiic major axes of velocity ellipsoids superimposed 
on the surface density distribution of the kinetic model K2 (Qs = 
1.3) at the time t = 1.3 Gyr. The angle between the major axis of 
velocity ellipsoid and the centre-anticentre direction {(f) = const) 
gives a value of vertex deviation. 

5 SUMMARY AND CONCLUSIONS 

We have developed the Boltzmann moment equation ap- 
proach for the dynamics of stars (BEADS-2D), which is a 
finite-difference Eulerian code based on the Boltzmann mo- 
ment equations up to second order. Wc adopt the zcro-licat- 
fiux approximation to close the system of moment equa- 
tions. The BEADS-2D code is formulated in the polar coor- 
dinates {r,4>), which makes it most suitable for the numeri- 
cal modelling of fiattened galactic stellar disks. By calculat- 
ing the tensor of the velocity dispersions, we can follow the 
anisotropy and the vertex deviation of the stellar system. 
Compared to numerical codes assuming a polytropic equa- 
tion of state (see e. g. Korchagin et al. 2000), the BEADS- 2D 
code allows for a larger variety of simulations, because the 
mass and velocity dispersion profiles can be chosen inde- 
pendently. For the reader's convenience, we provide the full 
set of three-dimensional Boltzmann moment equations up 
to second order in Appendix A . These equations are formu- 
lated in the cylindrical coordinates {z, r, (f>) using the usual 
zero-heat-flux approximation and are suitable for a global 
modelling of galactic disks without simplifying geometry as- 
sumptions. 

As an example of the utility of the BEADS-2D code, 
we study the time-dependent evolution of exponential stel- 
lar disks characterized initially by different values of the 
Toomre parameter Q^. The disks are embedded in a static 
dark matter halo, yielding a nearly fiat rotation curve (ex- 
cept for the innermost several kiloparsccs characterized by 
rigid rotation). Specifically, we find the following. 

(i) We confirm the results of linear stability analysis by 
Polyachenko et al. (1997), who obtained a stability thresh- 
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old value of Qs ~ 3.15 for stellar disks with a purely flat 
rotation curve. A physical interpretation for the instability 
is swing amplification of short-wavelength trailing distur- 
bances, which propagate through the disk centre and pro- 
vide a feedback mechanism for the swing amplifier. 

(ii) The character of the instability is distinct in the lin- 
ear and nonlinear regimes, the tentative boundary between 
which we define by the global Fourier amplitude C2{t) =0.1. 
As a reference model, we investigated a stellar disk with 
Qs = 1.3. In the early linear phase, instability demonstrates 
a typical lumpy pattern, which grows in amplitude with time 
and gradually transformes in to a global two-armed spiral 
in the late linear phase. This lumpy pattern results most 
probably from the interference between leading and trailing 
disturbances, propagating through the disk centre. A simi- 
lar phenomenon is discussed in Binney & Tremaine (1987) in 
connection with the numerical simulations of Toomre (1981). 
In the nonlinear phase, the redistribution of mass and an- 
gular momentum by the gravitational torques of spiral arms 
comes into play. The gravitational torques produce an inflow 
of mass, which forms a transient central bar inside corota- 
tion, and an outflow of mass, which forms a ring at the OLR. 
This process of mass redistribution contirmcs until the two- 
armed spiral is weakened by the shortage of matter. Finally, 
a transient bar is transformed into the compact dense cen- 
tral disk and the spiral axms virtually disappear, leaving 
only the outer diffuse ring at the OLR. 

(iii) The instability in disks with Qs ^ 2.5 is mostly 
dominated by the symmetric m = 2 mode, while the lop- 
sided m = 1 mode in disks with Qs > 2.5 may compete 
with or even prevail over the m — 2 mode. This is due 
the fact that the m = 1 mode has no inner Lindblad res- 
onance, which makes it difficult to stabilize. On the other 
hand, the stabilizing effect of the inner Lindblad resonance 
on the m = 2 mode is expected to grow along the sequence 
of increasing Qs. 

(iv) Stellar disks near a stability limit of Qs — 3.15 
still develop a clumpy structure. However, both the growth 
timescales (comparable to the Hubble time) and the small 
amplitudes of positive density perturbations make such in- 
stabilities difficult to observe at best. 

(v) Stellar disks are stabilized by one or more of the 
following procedures, a) An increase in the stellar velocity 
dispersions suppresses the growth of short-wavelength den- 
sity perturbations and thus inhibits the swing amplifier, b) 
Slower pattern speeds in hotter disk make it more likely for 
the inner Lindblad resonance to occur, thus providing the 
cutoff for the sing amplification feedback. 

The presence of the non-diagonal velocity dispersion 
terms in the Boltzmann moment equations allows us to 
study the vertex deviation. In most regions of the stellar disk 
the velocity ellipsoid is well aligned with the radial axes of 
the polar grid. In agreement with analytical predictions by 
Kuijkcn & Tremaine (1994), considerable vertex deviations 
show up in regions with strongly perturbed mass distribu- 
tions, i.e. near the spiral arms. The vertex deviations are 
large at the convex edge of the spiral arms, whereas they are 
small at the concave edge. The mean vertex deviations cor- 
relate well with the global Fourier amplitudes. Near the con- 
vex edge of the spiral arms, the ratio of radial to azimuthal 
components of the velocity ellipsoid can deviate considerably 
from the values predicted from the epicycle approximation. 



A more detailed numerical study of the properties of veloc- 
ity ellipsoids in spiral galaxies will be presented in a future 
paper. 
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APPENDIX A: THE 3D-BOLTZMANN 
MOMENT EQUATIONS UP TO SECOND 
ORDER 

The Boltzmann moment equations (up to second order) read 
in cylindrical coordinates (r, 0, z) as: 
Continuity equation: 

i + -a;:(-P-'-) + -^(p-.) + ^(p-«.) = o. (Ai) 

Momentum equations: 
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Velocity dispersion equations: 
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OfF-diagonal elements of the velocity-dispersion ten- 
sor: 
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Note that terms of third or higher order are neglected (zero- 
heat-flux approximation). 



APPENDIX B: THE SIX COMPONENTS OF 
THE ARTIFICIAL VISCOSITY STRESS 
TENSOR Q 

Here we give the components of the artificial viscosity stress 
tensor used in our simulations. 
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APPENDIX C: TESTS AND ACCURACY 

In this section we provide the results of two test problems 
probing the advection scheme and, especially, the conserva- 
tion of the specific angular momentum. The results of other 
tests will be given elsewhere. 



CI Test 1: relaxation problem 

The van Leer advection scheme was tested on a " relaxation" 
problem, in which a disk of constant density and velocity 
field proportional to r {ur = uor) is set, and the density is 
allowed to decrease. For this problem, the analytic solution 
of the continuity equation is E(t) — Soe^^""' and it gives 
E = 6.14 X 10^'' Eo at uot = 6. The numerical solution for 
this problem gives density values of 5.98 x 10~® Eo for the 
resolution of 256 radial grid points. When the density has 
decreased by nearly six orders of magnitude, its radial profile 
remains perfectly flat and the relative error is only 2.6%. 



C2 Test 2: Conservation of specific angular 
momentum 

An important point of concern when numerically studying 
the dynamics of galaxies is the ability of a code to conserve 
specific angular momentum. A comprehensive test problem 
on specific angular momentum conservation that covers ad- 
vection as well as pressure and gravity terms in the momen- 
tum equations was designed by Norman et al. (1980). These 
authors considered the isothermal gravitational collapse of 
rotating axisymmetric prestellar clouds in cylindrical geome- 
try. For a fluid with no mechanism for redistributing angular 
momentum, the mass M{K) in the cloud with specific angu- 
lar momentum K = Q,r^ less than or equal to if is a constant 
of motion. Deviations from the initial spectrum M{K) re- 
veal a redistribution of angular momentum. We consider a 
similar test problem of isothermal gravitational collapse of a 
flattened axisymmetric prestellar cloud (further referred as 
a "disk" ) that can be studied in the thin-disk approximation 
(see e.g. Basu 1997). This problem should test the ability of 
our polytropic model (for which we set 7 = 1 to switch to 
the isothermal regime) to conserve the speciflc angular mo- 
mentum. For a uniformly rotating gas disk (fi — Qq) with 
the radial surface density distribution given by Basu (1997), 



E(r) 



f'oEo 



the initial specific angular momentum spectrum is 



M{K) = 27rEorg 



'1 + 



K 



(CI) 



(02) 



where Kq = fio^'oi is the central surface density, and 
^0 = Cs/(1.5GEo) is the characteristic scale length which is 




K=£2r' 

Figure CI. The specific angular momentum spectrum of a col- 
lapsing flattened cloud. The quantity M{K) is the total mass 
in the cloud with specific angular momentum less than or equal 
to K = Clr^. The solid line shows the theoretical spectrum at 
t = yr. The filled circles gives M{K) at 0.185 Myr after the be- 
ginning of the simulation when the central protostar has already 
formed. A good agreement between both spectra indicates virtu- 
ally no angular momentum redistribution, as is indeed expected 
for the axisymmteric collapse. 



comparable to the Jeans length Aj = Cb/(GE) of a non- 
hoose Eo = 3. 55 
(equivalent to T 



rotating cloud. We choose Eo = 3.55 x 10 ^ gm cm ^, 



0.188 km s" 



10 K), fio 



s " pc" , and the mean molecular weight 2.33 
(molecular hydrogen with a 10% admixture of atomic he- 
lium). The disk has the inner and outer radii of nn = 10 AU 
and fout = 20000 AU. Wc introduce a "sink cell" at r < 
10 AU and impose a free inflow inner boundary condition 
and a constant mass and volume outer boundary condition. 
The mass of our disk is 2.45 M©. 

The ratios of rotational and thermal energies to the 
gravitational energy of the initial conflguration are 0.7% and 
25%, respectively. The disk is thus gravitationally unstable 
and is allowed to collapse under its own gravity. The initial 
theoretical spectrum M{K) is shown in Fig. CI by the solid 
line. The early evolution is characterised by a slow gravita- 
tional contraction and a subsequent runaway collapse. When 
the density in the sink cell exceeds 13.3 g cm^'^, the central 
protostar is assumed to form (due to a transition to the 
adiabatic phase). The flUed circles plot M{K) computed 
at t = 4300 yr after the formation of the central proto- 
star (0.185 Myr after the beginning of the simulation). At 
this time, approximately 0.3 M© has been accreted by the 
protostar. As is clearly seen, M{K) merges with the ini- 
tial spectrum (except for the end points where a minimal 
deviation is observed), indicating virtually no angular mo- 
mentum redistribution due to cither physical or rmmerical 
reasons. This is expected in the axisymmetric collapse with 
very little numerical diffusion. 



